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1  Introduction 


Numerical  evaluation  of  scattered  fields  generated  by  surfaces  of  revolution  requires 
fast  and  accurate  computation  of  individual  azimuthal  Fourier  components  Gm  of  the 
free-space  Green’s  function  for  the  Helmholtz  equation,  and  their  partial  derivatives 
(see,  for  example,  Harrington  and  Mautz  [9],  Medgyesi-Mitschang  and  Putnam  [11], 
Govind,  Wilton,  and  Glisson  [7],  Berton  and  Bills  [3],  Gedney  and  Mitra  [5],  Bon- 
nemason  and  Stupfel  [4],  Grannel,  Shirron,  and  Couchman  [8]).  The  functions  Gm 
are  defined  in  (10)  below. 

In  this  paper  we  derive  a  recurrence  relation  for  the  functions  Gm  which  (unlike 
three-term  recursion  formulae  for  most  classical  special  functions)  turns  out  to  involve 
five  terms.  In  addition  we  obtain  explicit  formulae  expressing  partial  derivatives  of 
these  functions  through  linear  combinations  of  GmS.  Such  properties  of  the  sequence 
[Gm]  are  analogous  to  that  of  classical  special  functions  and  orthogonal  polynomials 
(see,  for  example,  Abramovitz  and  Stegun  [1]),  and,  similarly  to  these  classical  objects, 
provide  for  fast  and  accurate  computation  of  an  individual  Gm  and  any  of  its  partial 
derivatives. 

The  plan  of  the  paper  is  as  follows.  In  Section  2  we  introduce  basic  notation  and 
definitions.  In  Section  3  we  obtain  the  generating  functions  for  the  sequence  [Gm] 
and  an  auxiliary  sequence  {Qm}  (see  its  definition  in  (13)  below).  In  Section  4  we 
derive  the  recurrence  relation  for  the  functions  Gm-  Iii  Section  5  we  express  partial 
derivatives  of  Gm  as  linear  combinations  of  these  functions.  In  Appendix  we  express 
certain  infinite  series  involving  functions  Gm  in  a  closed  form,  and  obtain  expansions 
of  these  functions  in  terms  of  spherical  Bessel  functions. 

2  Notation  and  Definitions 

The  free-space  Green’s  function  for  the  three  dimensional  Helmholtz  equation  ■0  sat¬ 
isfies  the  equation 

A0  +  A:V  = -y^(|r-r'l),  (1) 
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where  6  is  Dirac’s  delta  function,  r  and  r'  are  vectors  in 


k  27r/A, 


(2) 


and  A  is  the  wavelength. 

The  solution  of  (1)  satisfying  the  outgoing  radiation  condition  is  given  by  the 
formula 

exp(it|r-r-|)  ,3, 

(A-lr-r'l) 

Formulae  (l)-(3)  can  be  found  (in  a  slightly  different  form)  in  Chap.  7  of  Tikhonov 
and  Samarskii  [13]. 

In  polar  coordinates 


we  have 


where 


and 


r' = 

(4) 

^[r  —  r'p  =  a  —  bcos{<f>  —  4>'), 

(5) 

=  k^[p^+ipr+{z-zn 

(6) 

(7) 

z,  and  z'. 

a  >  b. 

(8) 

Combining  (3)  and  (5)  we  have 


_  ,,,  _  exp  (.-(g  -  i.cos(^  - 
nk\r  r|j-  • 


(9) 


The  functions  Gm  are  defines  as  Fourier  coefficients  of  (9)  with  respect  to  the 
azimuthal  angle  (j)  and  have  the  form 

exp  (i{a  —  b  cos 
{a  —  6  cos  4>y^'^ 


Gm{a,b)=  [  ^(A:|r-r'|)  •exp(fm^^)d(?!>=  /  ^  ^  ^  •  exp(?m(/))d<?!>. 

Jo  Jo 


(10) 


« 


I 
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Obviously,  for  any  a  >  b  and  integer  m, 


b)  —  G-mia,  i>). 


It  follows  from  (10)  that 


exp(i(a-6cos.^)'/2)  1 

/■  =  ^  E  G„(a,b)exp(rm^). 


/  ,  JM/2  0_  ^  - -ry  V*-/ 

(a  -  fecos(?i)i/2 

We  will  also  consider  the  sequence  of  functions  {Qm},  that  for  any  integer  m  and 
a  >  b  are  defined  by  the  formula 


(5m(a,6)=^  J  exp  —  6cos  •  exp(n 


which  implies  that 


1  ^ 

exp  (i(a  -  6cos  ^  Qm{a,b)exp{im(f>). 

m=— oo 


3  The  Generating  Functions 

The  main  result  of  this  section  is  Theorem  3.1  below  where  we  obtain  the  generating 
function  G  for  the  sequence  [Gm]  and  the  generating  function  Q  for  the  sequence 
It  is  well  known  that  generating  functions  are  normally  used  for  the  derivation 
of  recurrences  and  certain  other  properties  of  classical  special  functions  and  orthog¬ 
onal  polynomials  (see,  for  example,  Chap.  7  of  Lebedev  [10],  Chap.  11  of  Seaborn 
[12]).  The  corresponding  computations  for  the  sequence  [Gm]  require  analysis  of  two 
generating  functions  G  and  Q,  which  is  a  departure  from  the  standard  scheme  for 
classical  special  functions  where  usually  only  one  generating  function  is  involved. 

We  define  the  generating  function  P  for  a  certain  sequence  {Pm}  in  a  standard 
way,  i.  e.  as  a  function  whose  Laurent  expansion 

OO 

P{x,t)=  Y.  Pm{x)-t^  (15) 

m=— OO 

converges  uniformly  in  a  certain  ring  of  the  complex  t-plane. 


Theorem  3.1.  Let  for  any  a  >  b  the  functions  G  and  Q  be  defined  by  the  formulae 

.  ydef^  exp(^•(a-6(^-frl)/2)l/2) 

Gia.b.t)  =  2'K  •  - - ryr - ,  (lo) 
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and 


Q{a,  6,  t)  27r  •  exp  ^z(a  —  b{t  +  t  . 


(17) 


Then 

OO 

G{a,b,t)=  £  Gm{a,b)-t^,  (18) 

m=— OO 

and 

CO 

Q{a,b,t)=  ^  Q„(a,6)-r.  (19) 

m=— OO 

The  series  (18)  and  (19)  converge  uniformly  in  any  domain  lying  entirely  within  the 


ring 


t-  <  \t\<t+, 

(20) 

where 

(_  W  C-’  -  (c-=  -  l)  ,  4  =  c-'  +  (c-^  -  l)''’  . 

(21) 

and 

def  b 

c  =  -. 
a 

(22) 

1 

Proof.  Since  the  quadratic  equation 

• 

a-6(t  +  r')/2  =  0 

(23) 

has  the  roots  t_  and  it  is  easy  to  see  that  the  functions  G  and  Q  are  analytic 
in  the  ring  (20).  By  Laurent’s  theorem  each  of  them  has  a  Laurent  expansion  that 
converges  uniformly  in  any  domain  that  lies  wholly  in  this  ring.  Now  we  obtain  the 
formulae  (18)  and  (19)  by  substituting 

t  =  exp(t(;5)  (24) 


into  (12)  and  (14),  respectively.  • 

Remark  3.1.  Obviously,  the  functions  G  and  Q  are  connected  via  the  relation 

7  _ 0 _ 

which  follows  from  (16)  and  (17).  • 


(25) 
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4  The  Recurrence  Relation  for  the  Functions  Gm 


The  main  result  of  this  section  is  the  recurrence  relation  for  the  functions  Gm-  Namely, 
in  Theorem  4.1  below  we  show  that  for  all  a  >  6  and  |m|  >  1, 


aGm(O)  ~  2  “t"  ^))  ~  A*m(Gm+l(0)  b)  Gm-l{o-^  b))  + 


where 


def 

~  4m 


(26) 


(27) 


We  begin  with  two  preliminary  results  summarized  in  Lemma  4.1  and  Corollary  4.1 
below. 


Lemma  4.1.  For  any  integer  m  and  a  >  b, 

Q j  b^  —  i/J-m^Gm-i-l  Gm—l  (®)  ^))‘ 

Proof.  The  combination  of  (16)  and  (17)  yields 

.  (1  _  r^) .  G{a,  b,  t). 

Substituting  (18)  and  (19)  into  (29)  we  have 


(28) 


(29) 


oo  ,7,  so 

^  Q„{a,b)-m-r-'  =  -j  Y,  (G„(<i,6)-G„+2(o,6))-r.  (30) 

m  =  — OO  771=  — oo 

Equating  the  coefficients  of  equal  powers  of  t  in  (30)  we  immediately  obtain  (28).  • 
The  following  corollary  is  an  immediate  consequence  of  the  formula  (28). 

Corollary  4.1.  For  any  integer  |m|  >  1, 

^771+1(^7^)  Qm—li^ib') 


I  fir 


{ ^G„+7(<..6)  -  +  -!^G™-7(<..6))  .  .  (31) 

\m  +  l  m-' —  1  m  — 1  j 


Theorem  4.1.  For  any  a  >  b  and  m  >  1  the  relation  (26)  holds. 
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Proof.  We  start  with  the  identity 


1/2 


+  ^  )  •G(a,6,t)=  a--{t  +  t^) 


Q{a,bJ), 


(32) 


which  is  an  immediate  consequence  of  (25).  Differentiation  of  (32)  with  respect  to  t 
in  combination  with  (25)  yields 

dG{a,  b,  t) 


^)  '  G{a,b,t)  +  (a  — -{t  +  t  ^) 


dt 


-  J  ■  (1  -  r")  •  G(a,  b,t)-j-(l-r^)Q{a,b,t). 

Substituting  (18)  and  (19)  into  (33)  we  have 
b 

—  T  •  X/  {Gm{a,b)  —  Gm+2{<lib))  •  t"'  -\- 


(33) 


m  ■  b 


^  [a  ■  771  •  Gm{a,b) - —  •  (G^+i (a, b)  +  Gm-i (a, ft))  •  ^ 


m  — 1 


ib  ^ 

-T-  E  (Qn,(a,b)-Q„^,(a,b))-r.  (34) 

^  Tn=-oo 

Equating  the  coefficients  of  equal  powers  of  t  in  (34)  and  using  (31)  we  immediately 
obtain  (26).  • 

Remark  4.1.  An  important  for  applications  property  of  any  recurrence  relation  is 
its  numerical  stability,  i.  e.  the  sensitivity  of  the  {m  +  n)-th  term  of  the  sequence 
{m  =  0,±1,±2,...,  n  =  0,  ±1,  ±2, ...),  computed  via  n  steps  of  recursion,  to  small 
perturbations  of  its  m-th  term  (see,  for  example,  Chap.  9  of  Abramovitz  and  Stegun 
[1]).  A  somewhat  involved  analysis  shows  that  (26)  is  stable  for  both  upward  and 
downward  recurrences  if 

1  <  m  <  r_,  (35) 

is  stable  for  downward  and  unstable  for  upward  recurrences  if 

r_  <  m  <  r+,  (36) 

and  is  unstable  for  both  upward  and  downw'ard  recurrences  if 

m  >  r+.  (37) 


V 
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In  (35H37), 


1/2 


where  t-  and  t+  are  defined  in  (21).  The  proofs  of  these  results  will  be  reported 
elsewhere.  • 

Remark  4.2.  Simple  analysis  indicates  that  in  the  region  (37)  the  functions  Gm  are 
almost  always  (numerically)  small.  In  fact,  the  definition  (10)  and  Riemann  s  lemma 
yield 

lim  Gm(a)^)  =  0- 

m-+oo 

Combining  formulae  (9.3.1)  and  (10.1.1)  of  [1]  with  (67)  and  (68)  of  Appendix  it  is  easy 
to  show  that  the  decay  of  functions  Gm  for  m  satisfying  (37)  can  by  approximately 
described  by  the  formula 

/  \  m 

G„(<>,<')~(^)  ■ 

Therefore  the  functions  Gm  in  the  region  (37)  are  small  unless  r+  »  r_,  i.  e.  when 
a  K  b.  • 


(38) 


5  Partial  Derivatives  of  the  Functions  Gm 


In  this  section  we  show  that  partial  derivatives  of  functions  Gm  oan  be  expressed  via 
linear  combinations  of  these  functions. 


Lemma  5.1.  For  any  a  >  b, 


and 


dGm{a,b)  ^  2  .  p  .  (z  -  z')  • 

dz  ^  ’  da 


Furthermore,  for  any  integer  m, 


dGm{o.,b)  _  if  dGm+ijO’^b)  dGm-i{a,  b) 
m  “~2V  da  da 


(41) 


(42) 


(43) 


7 


Proof.  The  formulae  immediately  follow  from  (6)  and  (7).  The  formula  (43)  can  be 
easily  obtained  by  substituting  the  expansion  (18)  into  the  obvious  relation 


_1  dG{a,b,t)  dG{a,b,t) 

2'^^  ’'  da  db  ’ 


(44) 


and  equating  the  coefficients  of  equal  powers  of  t.  • 

Throughout  the  proofs  of  Lemma  5.2  and  Theorem  5.1  below  we  will  write  for 
brevity 


V 

■^m  — 


def  dG-mi^^^  b^ 


da 


(45) 


Lemma  5.2.  For  any  a  >  b  and  n  >  0, 


a\n  b\fi—i  —  Pm 


(46) 


where 


p„  t'  {n-l/2)-G„(a,b)-yp„- (G„^i(a,  b)  -  G„.,  (a,h)).  (47) 


Proof.  Substituting  the  expansion  (18)  into  the  obvious  relation 

b  ,,  dG{a,b,t)  dG{a,b,t) 

2'^  da  dt  ’ 

and  equating  the  coefficients  of  equal  powers  of  t  we  have 

-G'„(a,6)  =  pn  •  (-^n+l  ~  2f„_]). 

Next,  differentiation  of  (32)  with  respect  to  t  yields 


G(a.M)  + 


(48) 


(49) 


=  2^{a,b,l)+ -Q{a,b,t).  (50) 


Substituting  the  expansions  (18)  and  (19)  into  (50)  and  equating  the  coefficients  of 
equal  powers  of  t  we  have 


—  ^(^n+l  +  -^n-l)  +  ^  ^)  —  2  ’  ^)- 


(51) 


Combining  (28)  and  (51)  we  obtain 


-  ^G4a,b)  -  i//„(G„+i(a,6)  -  G„_i(a,6))  =  -  ^(X„+i  +  X„_a).  (52) 
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i 


4 


4 


« 


Now  (46)  is  an  immediate  consequence  of  (49)  and  (52).  • 

Theorem  5.1.  For  any  a  >  b  and  m  >2, 

dGm-2{a,b) _ An _ 

da  m{a?  —  —  6^  ’ 

where 

Am  *=  {2ma?  —  (m  —  1)5^)  •  {bpm-i  +  o,Pm)  ~ 

•  (m  —  1)  •  (5pm+i  +  flPm+2)  +  4  •  {rrY  —  1)  •  o'*  •  Sm,  (54) 


=  4(G™„(a,  6) -G„_, (<.,(.)),  (55) 

zb 

and  Pm  and  Pm  are  defined  in  (27)  and  (47),  respectively. 

Proof.  Differentiating  (26)  with  respect  to  a  we  obtain 

Gn{a,  b)  +  aXn  —  ^  (-^n+l  +  ^n-l)  “  Pn{Xn+l  “  Xn-l)  + 


”  Xn+2--4^^n  +  -^A„_2)  =0. 


+  n2_i-—  n-1 

Substituting  (49)  and  (52)  into  (56)  we  have 

1 


n  —  1 


A„_2- 


2^  Y  I  ^  Y  —  " 

Zi-n  T - rT^n+2  — 


n 


2  -  r 


n  +  1 


(56) 


(57) 


where  s„  is  defined  in  (55). 

The  combination  of  (46)  forn  =  m  -  1,  m,  m  +  1  and  m  +  2  with  (57)  forn  =  m 
yields  the  following  system  of  linear  equations: 

aXm-l  —  bXm—2  —  Pm-\i 
aXm  bXm—1  —  Prm 


CtAjn+1  bXm  —  Pm+l5 
aXm-{-2  bXm+1  Pto+2? 
1 


‘Xm  —  2 


2^  Y  I  ^  Y  _  o 

■^m  - i  T-^‘m+2  — 


-  -iTn-2  “  9  ^  I  1 

m  —  1  —  1  m  +  1 


(58) 
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Evaluating  ^”,71-2  from  (58)  and  using  the  definition  (45)  we  obtain  (53)“(55).  • 

Remark  5.1.  Clearly,  for  all  integer  m  and  a  >  b  any  partial  derivative  d^^^Gmldp^dz^ 
(it,  I  =  0, 1,2...)  can  be  expressed  as  a  linear  combination  of  a  finite  number  of  func¬ 
tions  Gm,  which  is  a  consequence  of  (11),  (41)-(43),  and  (53)-(55).  • 

6  Conclusions 

We  have  demonstrated  that  the  functions  Gm  satisfy  a  recurrence  relation  and  showed 
that  their  partial  derivatives  can  be  expressed  via  linear  combinations  of  these  func¬ 
tions  themselves.  Such  properties  of  the  sequence  of  functions  {Gm}  ^^re  similar  to 
that  of  all  the  families  of  classical  special  functions  and  orthogonal  polynomials. 

The  use  of  the  recurrence  (26)  significantly  reduces  the  complexity  of  the  numer¬ 
ical  computation  of  the  scattered  fields  generated  by  surfaces  of  revolution.  As  was 
observed  by  many  authors  such  computations  can  be  reduced  to  the  solution  of  M 
integral  equations  on  the  generating  curve  of  the  scatterer  (see,  for  example,  [7]). 
Here  M  ~  A:  is  the  number  of  nonvanishing  (with  the  given  accuracy)  azimuthal 
Fourier  components  of  the  incident  field  computed  on  the  surface  of  the  scatterer, 
and  k  is  defined  in  (2).  After  appropriate  discretization  these  integral  equations  are 
reduced  to  a  sequence  of  linear  systems  with  dense  matrices  (m  =  0, 1,  ...,M). 
Any  element  of  these  matrices  is  a  linear  combination  of  a  finite  number  of  the 
functions  Gm  and  their  partial  derivatives,  with  the  arguments  of  these  functions  p, 
2,  p',  and  z'  belonging  to  the  mesh  on  the  generating  curve  of  the  scatterer  (see,  for 
example,  [3],  [5]). 

It  is  well  known  that  evaluating  the  matrix  elements  consumes  a  major 

portion  of  the  CPU  time  (see,  for  example,  [5]).  Indeed,  simple  arguments  based  on 
the  Nyquist  theorem  show  that  the  number  of  elements  in  Z^'”)  is  at  least  0{N'^), 
where  N  is  the  size  of  the  scatterer  in  wavelengths.  Note  that  N  ^  k,  which  is  a 
consequence  of  (2).  The  computation  of  one  Z,-J”^  requires  0{N)  operations  which 
immediately  follows  from  (10).  Therefore  the  CPU  time  Ti  for  computing  all  the 
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matrix  elements  has  the  estimate 


Ti  =  0  (m  •  N^)  .  (59) 

Next,  solving  the  linear  system  for  every  m  by  means  of  some  iterative  technique 
normally  requires  application  of  the  system’s  matrix  to  a  certain  sequence  of  vectors 
(see,  for  example.  Chap.  10  of  Golub  and  Van  Loan  [6]),  i.  e.  can  be  done  for  0{N^) 
operations.  Thus  the  CPU  time  Tsoi  for  solving  the  M  linear  systems  has  the  estimate 

Tsoi  =  0{M-N^).  (60) 

Comparing  (59)  and  (60)  we  see  that  the  complexity  of  computing  the  scattered  field 
is  dominated  by  the  cost  of  generating  matrices  (m  ==  0, 1,  ...,M)  and  thus  has 
the  CPU  time  estimate  Ti. 

Now  we  turn  to  estimating  the  CPU  time  for  the  generation  of  the  matrix  elements 
by  means  of  (26);  this  method  consists  of  two  stages.  On  the  first  stage  we  must 
evaluate  the  starting  values  of  the  recursion,  i.  e.  0  (V^)  elements  of  the  four  initial 
matrices,  which  requires  0  (N^)  operations.  On  the  second  stage  we  compute  all  other 
matrices  using  the  recursion  (26)  which  requires  0  {M  •  N'^)  operations.  Therefore 
the  CPU  time  T2  of  computing  the  matrices  by  means  of  (26)  has  the  estimate 

T2  =  0  [N^)  +0{M-N^).  (61) 

Obviously,  now  the  CPU  time  estimate  of  the  computation  of  the  scattered  field  has 
the  form  (61)  which  immediately  follows  from  (60)  and  (61). 

For  large-scale  problems  we  have  iV  >  1  and  M  >  1,  which  in  combination  with 
(59)  and  (61)  yields 

T2  <  Ti.  (62) 

Currently  a  fast  algorithm  for  generating  matrices  using  (26)  is  being  imple¬ 
mented  and  the  corresponding  results  will  be  reported  at  a  later  date. 
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8  Appendix 


In  Appendix  we  discuss  certain  consequences  of  the  formulae  (10)  and  (12). 
Theorem  A.l.  For  any  a  >  b, 

exp  —  by^^'^ 


Gm(o,  b)  —  Go{a,  b)  +  2  Gmjo-,  b)  —  217- 


m—\ 


{-\rGm{a,  b)  =  Goia,  b)  +  2Y  i-irGmia,  b)  =  2x 

a  =  -oo  m=l 

oo  oo 

Yi  |G„(o,6)P=  |G„(a.6)p  +  2  5:  |G„(a,6)]^  = 


[a  -  by/'^  ’ 

exp  (y(a  +  by^^'j 


{a  +  6)1/2 
47r^ 


(63) 


(64) 


(65) 


fa2_i2')l/2- 

m  =  —  oo  m=l  ^  ’ 

Proof.  Formulae  (63)  and  (64)  immediately  follow  from  (12)  after  substituting  (f)  =  0 
and  ^  —  n,  respectively.  The  formula  (65)  is  easily  obtained  by  applying  Parseval’s 
theorem  to  (12)  and  using  the  relation 

'■2’^  d<j)  277 


I- 

Jo  a 


(66) 


bcos(f)  —  b^y/'^ 

Theorem  A. 2.  For  any  a  >  b  and  integer  I  >  0, 

OO 

G2t{a,  b)  =  2i  X](47i  +  1)  •  A(7?  -  /)  •  A(n  +  /)  •  j2nir-)  ■  h[^^{r+),  (67) 


n=/ 


G2l+l{a,b)  =  27^^(477  +  3)  •  A(77  -  1)  •  A{n  +  /  +  1)  ■  j2n+](r-)  •  62^„+i(7’+),  (68) 


n-=l 


where  and  j  and  are  spherical  Bessel  functions  of  the  first  and  third  kind,  respec¬ 
tively,  r^  and  r_  are  given  in  (38),  and  A  is  defined  by  the  formula 

Ww  <••> 

Proof.  Combining  formulae  (10.1.1),  (10.1.45),  and  (10.1.46)  of  [1]  we  have 
exp  (zYa  —  6cos  <6)1/^ )  ^ 

- (a  _  6 cos  <l>y/^ - ^  +  1)  •  i„(r_)  •  hy\r+)  ■  P„(cos  fi),  (70) 

where  Pn  are  Legendre  polynomials.  Following  Alpert  and  Rokhlin  [2]  we  can  write 
for  any  integer  p  >  0 


1 


P2p(cos  <f>)  =  -  ^(2  -  6og)  •  A(p  -  q)  ■  A(p  +  q)  ■  cos{2q(f>), 


77 


(71) 


g=0 


12 


* 


4 


4 


and 

-p2p+i(cos  <?!>)  =  —  51  ~  '  ^^P  "t"  9  +  1)  ■  cos((2g  +  1)<^),  (72) 

g=0 

where  Snm  is  Kronecker’s  delta.  Now  substituting  (70),  (71),  and  (72)  into  (10)  we 
immediately  obtain  (67)  and  (68).  • 

The  formula  (67)  for  /  =  0  is  given  in  [3]. 


Remark  A.l.  Obviously,  the  function  A  is  closely  related  to  the  beta  function  B 
(see,  for  example.  Chap.  6  of  [1]).  In  fact,  comparing  the  formula  (6.2.2)  of  [1]  with 
(69)  we  have 


K{x) 


B{x,ll2) 

r(i/2)  • 


Note,  that  the  function  A  satisfies  the  recurrence  relation 


A(0)  =  r(l/2)  =  7r^/^  A(n  +  l)  =  A(n)- 


n  +  1/2 
n  +  1 


for  all  n  =  0, 1, 


(73) 


(74) 


which  is  an  immediate  consequence  of  (69).  • 
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